knitr::opts_chunk$set(
echo = TRUE,
message = FALSE,
warning = FALSE)
library(tidyverse)
library(skimr)
library(lubridate)
library(stringr)
library(nortest)
library(moments)
library(rstatix)
library(effectsize)
library(igraph)
library(tidygraph)
library(ggraph)
library(ggiraph)
library(broom)
library(robustbase)
library(lmtest)
library(performance)
library(see)
library(modelsummary)
# Formázott p-érték függvény
format_p <- function(p) {
ifelse(p < 0.001, "< 0.001", as.character(round(p, 3)))}
A gyógyszerek forgalomba hozatali engedélyezése összetett, több lépcsős szabályozási folyamat, amelyet az engedélyezést követően is folyamatos felügyelet és a dokumentáció ismtételt aktualizálása kísér. A forgalomba hozatali engedélyhez tartozó dokumentáció módosításai, az úgynevezett revíziók, a gyógyszer teljes életciklusa során tükrözhetik többek között a biztonsági jelzések, a klinikai bizonyítékok bővülésének, illetve a szabályozói eljárások sajátosságainak hatását. Emellett az engedélyezések mennyisége és összetétele (például terápiás területek szerint) időben is változhat, ami az gyógyszerfejlesztési trendek és a gyógyszercégek fókuszának változását mutathatják.
Az Európai Gyógyszerügynökség (European Medicines Agency, EMA) központi szerepet játszik a gyógyszerek engedélyezésében az Európai Unióban. Jelen elemzés célja az EMA gyógyszerengedélyezési adatbázisának részletes vizsgálata, amely négy fő kutatási kérdés köré szerveződik: 1. Az engedélyezések időbeli trendjeinek feltárása 1995 és 2023 között. 2. Az ATC-csoportok (Anatomical Therapeutic Chemical) közötti különbségek vizsgálata kulcsváltozók mentén 3. A revíziók számát befolyásoló tényezők azonosítása többváltozós modellezéssel 4. A központi idegrendszeri (ATC: N csoport) készítményeknél a terápiás területek és hatóanyagok kapcsolati szerkezetének hálózati szemléletű bemutatása.
Az elemzéshez az EMA nyilvánosan elérhető adatbázisát használtam, amely a TidyTuesday projekt keretében került publikálásra (2023. március 14.).
Az elemzés reprodukálható módon, webes forrásból betöltött adatbázison alapul, közvetlenül a TidyTuesday GitHub tárhelyéről került beolvasásra.
# Adat beolvasása a GitHub-ról
drugs <-
readr::read_csv('https://raw.githubusercontent.com/rfordatascience/tidytuesday/main/data/2023/2023-03-14/drugs.csv')
# Adatok áttekintése
dim(drugs)
## [1] 1988 28
names(drugs)
## [1] "category"
## [2] "medicine_name"
## [3] "therapeutic_area"
## [4] "common_name"
## [5] "active_substance"
## [6] "product_number"
## [7] "patient_safety"
## [8] "authorisation_status"
## [9] "atc_code"
## [10] "additional_monitoring"
## [11] "generic"
## [12] "biosimilar"
## [13] "conditional_approval"
## [14] "exceptional_circumstances"
## [15] "accelerated_assessment"
## [16] "orphan_medicine"
## [17] "marketing_authorisation_date"
## [18] "date_of_refusal_of_marketing_authorisation"
## [19] "marketing_authorisation_holder_company_name"
## [20] "pharmacotherapeutic_group"
## [21] "date_of_opinion"
## [22] "decision_date"
## [23] "revision_number"
## [24] "condition_indication"
## [25] "species"
## [26] "first_published"
## [27] "revision_date"
## [28] "url"
glimpse(drugs)
## Rows: 1,988
## Columns: 28
## $ category <chr> "human", "human", "human",…
## $ medicine_name <chr> "Adcetris", "Nityr", "Ebva…
## $ therapeutic_area <chr> "Lymphoma, Non-Hodgkin; H…
## $ common_name <chr> "brentuximab vedotin", "ni…
## $ active_substance <chr> "brentuximab vedotin", "ni…
## $ product_number <chr> "002455", "004582", "00457…
## $ patient_safety <lgl> FALSE, FALSE, FALSE, FALSE…
## $ authorisation_status <chr> "authorised", "authorised"…
## $ atc_code <chr> "L01XC12", "A16AX04", NA, …
## $ additional_monitoring <lgl> FALSE, FALSE, TRUE, TRUE, …
## $ generic <lgl> FALSE, TRUE, FALSE, FALSE,…
## $ biosimilar <lgl> FALSE, FALSE, FALSE, FALSE…
## $ conditional_approval <lgl> FALSE, FALSE, FALSE, FALSE…
## $ exceptional_circumstances <lgl> FALSE, FALSE, TRUE, FALSE,…
## $ accelerated_assessment <lgl> FALSE, FALSE, FALSE, FALSE…
## $ orphan_medicine <lgl> TRUE, FALSE, TRUE, FALSE, …
## $ marketing_authorisation_date <date> 2012-10-25, 2018-07-26, 2…
## $ date_of_refusal_of_marketing_authorisation <date> NA, NA, NA, NA, NA, NA, N…
## $ marketing_authorisation_holder_company_name <chr> "Takeda Pharma A/S", "Cycl…
## $ pharmacotherapeutic_group <chr> "Antineoplastic agents", "…
## $ date_of_opinion <date> 2012-07-19, 2018-05-31, 2…
## $ decision_date <date> 2022-11-17, 2023-03-10, 2…
## $ revision_number <dbl> 34, 4, 2, 3, 30, 24, 4, 18…
## $ condition_indication <chr> "Hodgkin lymphomaAdcetris …
## $ species <chr> NA, NA, NA, NA, NA, NA, NA…
## $ first_published <dttm> 2018-07-25 13:58:00, 2018…
## $ revision_date <dttm> 2023-03-13 11:52:00, 2023…
## $ url <chr> "https://www.ema.europa.eu…
skimr::skim(drugs)
| Name | drugs |
| Number of rows | 1988 |
| Number of columns | 28 |
| _______________________ | |
| Column type frequency: | |
| character | 13 |
| Date | 4 |
| logical | 8 |
| numeric | 1 |
| POSIXct | 2 |
| ________________________ | |
| Group variables | None |
Variable type: character
| skim_variable | n_missing | complete_rate | min | max | empty | n_unique | whitespace |
|---|---|---|---|---|---|---|---|
| category | 0 | 1.00 | 5 | 10 | 0 | 2 | 0 |
| medicine_name | 0 | 1.00 | 3 | 125 | 0 | 1976 | 0 |
| therapeutic_area | 285 | 0.86 | 4 | 400 | 0 | 669 | 0 |
| common_name | 4 | 1.00 | 4 | 220 | 0 | 1261 | 0 |
| active_substance | 1 | 1.00 | 4 | 823 | 0 | 1345 | 0 |
| product_number | 0 | 1.00 | 6 | 6 | 0 | 1932 | 0 |
| authorisation_status | 1 | 1.00 | 7 | 10 | 0 | 3 | 0 |
| atc_code | 28 | 0.99 | 3 | 18 | 0 | 1074 | 0 |
| marketing_authorisation_holder_company_name | 4 | 1.00 | 4 | 65 | 0 | 615 | 0 |
| pharmacotherapeutic_group | 34 | 0.98 | 7 | 174 | 0 | 365 | 0 |
| condition_indication | 12 | 0.99 | 18 | 7597 | 0 | 1886 | 0 |
| species | 1709 | 0.14 | 4 | 67 | 0 | 59 | 0 |
| url | 0 | 1.00 | 53 | 148 | 0 | 1988 | 0 |
Variable type: Date
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| marketing_authorisation_date | 60 | 0.97 | 1995-10-20 | 2023-02-20 | 2013-06-09 | 1127 |
| date_of_refusal_of_marketing_authorisation | 1913 | 0.04 | 2004-09-07 | 2022-04-29 | 2013-04-25 | 67 |
| date_of_opinion | 779 | 0.61 | 1995-07-12 | 2022-12-15 | 2016-07-21 | 389 |
| decision_date | 45 | 0.98 | 1998-08-20 | 2023-03-10 | 2022-02-16 | 815 |
Variable type: logical
| skim_variable | n_missing | complete_rate | mean | count |
|---|---|---|---|---|
| patient_safety | 0 | 1 | 0.01 | FAL: 1977, TRU: 11 |
| additional_monitoring | 0 | 1 | 0.19 | FAL: 1601, TRU: 387 |
| generic | 0 | 1 | 0.16 | FAL: 1673, TRU: 315 |
| biosimilar | 0 | 1 | 0.05 | FAL: 1896, TRU: 92 |
| conditional_approval | 0 | 1 | 0.02 | FAL: 1940, TRU: 48 |
| exceptional_circumstances | 0 | 1 | 0.02 | FAL: 1940, TRU: 48 |
| accelerated_assessment | 0 | 1 | 0.02 | FAL: 1940, TRU: 48 |
| orphan_medicine | 0 | 1 | 0.08 | FAL: 1826, TRU: 162 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| revision_number | 96 | 0.95 | 13.53 | 11.65 | 0 | 4.75 | 11 | 19 | 89 | ▇▃▁▁▁ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| first_published | 0 | 1.00 | 1998-08-20 00:00:00 | 2023-03-09 18:50:00 | 2018-04-14 23:14:30 | 1760 |
| revision_date | 29 | 0.99 | 2000-07-17 02:00:00 | 2023-03-13 11:52:00 | 2022-05-11 11:58:00 | 1932 |
# Folytonos változók kialakítása
## Gyógyszer piacon töltött idejének (days_on_market változó) kiszámítása napokban (Az adatbázis feltöltési dátuma: 2023.03.13.)
drugs <- drugs %>%
mutate(
days_on_market = as.numeric(ymd("2023-03-13") - ymd(marketing_authorisation_date)))
## Ellenörzés - days_on_market
summary(drugs$days_on_market)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 21 1818 3564 3948 5655 10006 60
## Eloszlás - days_on_market
ggplot(drugs, aes(x = days_on_market)) +
geom_histogram(bins = 30, fill = "steelblue", color = "white") +
labs(x = "Napok száma a piacon", y = "Gyakoriság") +
theme_minimal()
A „piacon töltött időt” (days_on_market változó) a forgalomba hozatali engedély dátuma és az adatbázis feltöltési dátuma (2023-03-13) közötti napok száma alapján képeztem, ezzel biztosítva, hogy a mutató ne függjön a futtatás aktuális dátumától.
Az adatbázis összesen 1988 gyógyszert tartalmaz, amelyek közül 1706 (85.8%) humán és 282 (14.2%) állatgyógyászati készítmény. Az adatbázis 28 változót tartalmaz, beleértve a gyógyszer nevét, terápiás területét, hatóanyagát, ATC kódját, engedélyezési dátumát és státuszát, valamint olyan különböző jellemzőket a gyógyszerekkel kapcsolatban, mint pl. generikus, orphan, gyorsított eljárás stb.
A gyógyszerek engedélyezési dátuma 1995. október 20. és 2023. február 20. között oszlik meg (medián: 2013. június 9.). Az engedélyezési státusz tekintetében a gyógyszerek 79.1%-a (n = 1573) rendelkezik érvényes engedéllyel, 18.0% (n = 357) visszavont, és 2.9% (n = 57) elutasított státuszú.
A speciális státuszok tekintetében a gyógyszerek 15.8%-a (n = 315) generikus, 8.2%-a (n = 162) orphan gyógyszer és 2.4%-a (n = 48) részesült gyorsított eljárásban.
A gyógyszerek piacon töltött idejének eloszlása a teljes mintában széles tartományt fed le (Min: 21 nap; Mdn: 3564 nap; M: 3948 nap; Max: 10006 nap; NA: 60).
# Kategórikus változók gyakorisága és százalékos eloszlása
vars <- c(
"category", "patient_safety", "authorisation_status", "generic",
"conditional_approval", "exceptional_circumstances",
"accelerated_assessment", "orphan_medicine")
cat_tables <- drugs %>%
select(all_of(vars)) %>%
lapply(function(x) {
tab <- table(x, useNA = "ifany")
pct <- prop.table(tab) * 100
data.frame(
value = names(tab),
n = as.integer(tab),
percent = round(as.numeric(pct), 1),
row.names = NULL)})
cat_tables$category
## value n percent
## 1 human 1706 85.8
## 2 veterinary 282 14.2
cat_tables$patient_safety
## value n percent
## 1 FALSE 1977 99.4
## 2 TRUE 11 0.6
cat_tables$authorisation_status
## value n percent
## 1 authorised 1573 79.1
## 2 refused 57 2.9
## 3 withdrawn 357 18.0
## 4 <NA> 1 0.1
cat_tables$generic
## value n percent
## 1 FALSE 1673 84.2
## 2 TRUE 315 15.8
cat_tables$conditional_approval
## value n percent
## 1 FALSE 1940 97.6
## 2 TRUE 48 2.4
cat_tables$exceptional_circumstances
## value n percent
## 1 FALSE 1940 97.6
## 2 TRUE 48 2.4
cat_tables$accelerated_assessment
## value n percent
## 1 FALSE 1940 97.6
## 2 TRUE 48 2.4
cat_tables$orphan_medicine
## value n percent
## 1 FALSE 1826 91.9
## 2 TRUE 162 8.1
# Piacon töltött idő (days_on_market) leíró statisztikája
summary(drugs$days_on_market)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 21 1818 3564 3948 5655 10006 60
# Gyógyszerek engedélyezésének száma évente
drugs %>%
filter(!is.na(marketing_authorisation_date)) %>%
mutate(year = year(marketing_authorisation_date)) %>%
count(year) %>%
ggplot(aes(year, n)) +
geom_line() +
scale_x_continuous(breaks = seq(1995, 2022, by = 1), minor_breaks = NULL, limits = c(NA, 2022)) +
labs(title = "Gyógyszer-engedélyezések száma évente", x = "Év", y = "n") +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Top 10 gyártó
drugs %>%
count(marketing_authorisation_holder_company_name, sort = TRUE) %>%
head(10)
## # A tibble: 10 × 2
## marketing_authorisation_holder_company_name n
## <chr> <int>
## 1 Accord Healthcare S.L.U. 58
## 2 Novartis Europharm Limited 58
## 3 Pfizer Europe MA EEIG 43
## 4 Zoetis Belgium SA 40
## 5 AstraZeneca AB 36
## 6 Boehringer Ingelheim Vetmedica GmbH 36
## 7 Merck Sharp & Dohme B.V. 32
## 8 Teva B.V. 32
## 9 Intervet International BV 31
## 10 Eli Lilly Nederland B.V. 30
Az 1. ábra a gyógyszer-engedélyezések számának alakulását mutatja az 1995-2022 közötti időszakban. A trendvonal jelentős fluktuációt mutat, kiugró csúcsokkal 2001 (n ≈ 90), 2017 (n ≈ 110) és 2021 körül (n ≈ 120). A COVID-19 pandémia időszakában (2020-2021) a gyógyszer-engedélyezések száma jelentősen megugrott, ami feltételezhetően a vakcinák és vírusellenes szerek gyorsított engedélyezésének tudható be elsősorban.
A forgalmazási engedély jogosultjai közül a legtöbb gyógyszert az Accord Healthcare S.L.U. és a Novartis Europharm Limited jegyzi (egyaránt n = 58), amelyeket a Pfizer Europe MA EEIG (n = 43) és a Zoetis Belgium SA (n = 40) követ.
A humán készítmények ATC-főcsoport szerinti bontása lehetővé teszi annak áttekintését, hogy az engedélyezések elsősorban mely terápiás területekhez kapcsolódtak, illetve hogyan alakultak ezen időszakban. Az ATC-főcsoportok szerinti évekre vontott ábra a főcsoportok közötti nagyságrendi különbségeket és az időbeli dinamikát együttesen teszi láthatóvá, így a későbbi, célzottabb összehasonlítások megalapozására szolgálhat.
# Human gyógyszerek ATC kód csoportonként
drugs %>%
filter(category == "human", !is.na(marketing_authorisation_date), !is.na(atc_code)) %>%
mutate(
year = year(marketing_authorisation_date),
atc_group = str_extract(atc_code, "^[A-Z]")) %>%
count(year, atc_group) %>%
ggplot(aes(year, n)) +
geom_line() +
facet_wrap(~atc_group, scales = "free_x") +
scale_x_continuous(breaks = seq(1995, 2022, by = 5), limits = c(1995, 2022), minor_breaks = NULL) +
labs(
title = "Humán gyógyszer-engedélyezések száma évente ATC-rendszer szerint",
x = "Év",
y = "n",
caption = "Note. ATC besorolás: A - Tápcsatorna és anyagcsere; B - Vér és vérképzőszervek; C - Cardiovascularis rendszer; D - Bőrgyógyászati készítmények; G - Urogenitalis rendszer és nemi hormonok;\nH - Szisztémás hormonkészítmények; J - Szisztémás fertőzésellenes szerek; L - Daganatellenes és immunmoduláns szerek; M - Váz- és izomrendszer; N - Központi idegrendszer;\nP - Parazitaellenes készítmények; R - Légzőrendszer; S - Érzékszervek; V - Egyéb") +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
strip.text = element_text(size = 10),
plot.caption = element_text(hjust = 0, size = 9))
# Heatmap - ATC csoportok évenete
drugs %>%
filter(category == "human", !is.na(marketing_authorisation_date), !is.na(atc_code)) %>%
mutate(
year = year(marketing_authorisation_date),
atc_group = str_extract(atc_code, "^[A-Z]")) %>%
filter(year <= 2022) %>%
count(year, atc_group) %>%
group_by(year) %>%
mutate(percent = n / sum(n) * 100) %>%
ungroup() %>%
ggplot(aes(x = atc_group, y = factor(year), fill = percent)) +
geom_tile(color = "white") +
geom_text(aes(label = ifelse(n > 0, n, "")), size = 2.5) +
scale_fill_gradient(low = "white", high = "steelblue", na.value = "white") +
labs(
title = "Humán gyógyszer-engedélyezésének száma ATC-rendszer szerint",
x = "ATC csoport",
y = "Év",
fill = "Arány (%)",
caption = "Note. ATC besorolás: A - Tápcsatorna és anyagcsere; B - Vér és vérképzőszervek; C - Cardiovascularis rendszer; D - Bőrgyógyászati készítmények; G - Urogenitalis rendszer és nemi hormonok; \nH - Szisztémás hormonkészítmények; J - Szisztémás fertőzésellenes szerek; L - Daganatellenes és immunmoduláns szerek; M - Váz- és izomrendszer; N - Központi idegrendszer; \nP - Parazitaellenes készítmények; R - Légzőrendszer; S - Érzékszervek; V - Egyéb") +
theme_minimal() +
theme(
plot.caption = element_text(hjust = 0, size = 9),
panel.grid = element_blank())
# ATC csoportok gyakorisága a teljes időszakban
drugs %>%
filter(category == "human", !is.na(atc_code)) %>%
mutate(atc_group = str_extract(atc_code, "^[A-Z]")) %>%
count(atc_group, sort = TRUE) %>%
mutate(percent = round(n / sum(n) * 100, 1))
## # A tibble: 14 × 3
## atc_group n percent
## <chr> <int> <dbl>
## 1 L 463 27.5
## 2 J 253 15.1
## 3 A 201 12
## 4 N 186 11.1
## 5 B 135 8
## 6 C 95 5.7
## 7 V 72 4.3
## 8 R 69 4.1
## 9 G 59 3.5
## 10 M 54 3.2
## 11 H 40 2.4
## 12 S 34 2
## 13 D 18 1.1
## 14 P 2 0.1
Az ATC-rendszer szerinti elemzés a humán gyógyszerekre korlátozódott. A 14 fő ATC csoportból a leggyakoribb az L csoport (daganatellenes és immunmoduláns szerek; n = 463, 27.5%), amelyet a J csoport (szisztémás fertőzésellenes szerek; n = 253, 15.1%), az A csoport (tápcsatorna és anyagcsere; n = 201, 12.0%) és az N csoport (központi idegrendszer; n = 186, 11.1%) követ. Ez a négy csoport együttesen a humán gyógyszerek 65.7%-át teszi ki.
A heatmap-vizualizáció (2. ábra) az ATC csoportok éves megoszlását szemlélteti. Az L csoport dominanciája különösen az elmúlt évtizedben vált hangsúlyossá, ami feltehetően a precíziós onkológia és az immunterápiák fejlődésével magyarázható.
A további elemzéseket a négy leggyakoribb ATC csoportra (L, J, A, N) szűkítettem, amelyek együttesen 1103 gyógyszert foglalnak magukban.
# Top 4 ATC csoport kiszűrése
top4_atc <- drugs %>%
filter(category == "human", !is.na(atc_code)) %>%
mutate(atc_group = str_extract(atc_code, "^[A-Z]")) %>%
filter(atc_group %in% c("L", "J", "A", "N"))
# Leíró statisztika - revision_number
top4_atc %>%
group_by(atc_group) %>%
summarise(
n = n(),
hianyzo = sum(is.na(revision_number)),
atlag = mean(revision_number, na.rm = TRUE),
sd = sd(revision_number, na.rm = TRUE),
min = min(revision_number, na.rm = TRUE),
q25 = quantile(revision_number, 0.25, na.rm = TRUE),
median = median(revision_number, na.rm = TRUE),
q75 = quantile(revision_number, 0.75, na.rm = TRUE),
max = max(revision_number, na.rm = TRUE))
## # A tibble: 4 × 10
## atc_group n hianyzo atlag sd min q25 median q75 max
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 A 201 7 14.1 9.51 0 6 14 21 37
## 2 J 253 6 17.7 14.2 0 7 14 26.5 68
## 3 L 463 31 14.5 12.9 0 5 11 20 89
## 4 N 186 8 15.0 11.6 0 7 13 21 58
# Leíró statisztika - days_on_market
top4_atc %>%
group_by(atc_group) %>%
summarise(
n = n(),
hianyzo = sum(is.na(days_on_market)),
atlag = mean(days_on_market, na.rm = TRUE),
sd = sd(days_on_market, na.rm = TRUE),
min = min(days_on_market, na.rm = TRUE),
q25 = quantile(days_on_market, 0.25, na.rm = TRUE),
median = median(days_on_market, na.rm = TRUE),
q75 = quantile(days_on_market, 0.75, na.rm = TRUE),
max = max(days_on_market, na.rm = TRUE))
## # A tibble: 4 × 10
## atc_group n hianyzo atlag sd min q25 median q75 max
## <chr> <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 A 201 4 4182. 2524. 94 2163 3979 5908 9813
## 2 J 253 2 4300. 2774. 98 2008. 4042 6208. 9805
## 3 L 463 20 3130. 2455. 26 1138 2482 4676 9968
## 4 N 186 8 3987. 2420. 126 2267 3808. 5258. 9772
# ATC csoport nevek hozzáadása
top4_atc <- top4_atc %>%
mutate(atc_name = case_when(
atc_group == "L" ~ "L - Daganatellenes és\nimmunomoduláns szerek",
atc_group == "J" ~ "J - Szisztémás\nfertőzésellenes szerek",
atc_group == "A" ~ "A - Tápcsatorna\nés anyagcsere",
atc_group == "N" ~ "N - Központi\nidegrendszer"
))
# Violin plot - revision_number
ggplot(top4_atc, aes(x = atc_name, y = revision_number, fill = atc_name)) +
geom_violin(alpha = 0.5) +
geom_boxplot(width = 0.2, outlier.shape = NA) +
stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "red") +
labs(
title = "Revíziók száma ATC csoportonként",
x = "ATC csoport",
y = "Revíziók száma"
) +
theme_minimal() +
theme(legend.position = "none")
# Violin plot - days_on_market
ggplot(top4_atc, aes(x = atc_name, y = days_on_market, fill = atc_name)) +
geom_violin(alpha = 0.5) +
geom_boxplot(width = 0.2, outlier.shape = NA) +
stat_summary(fun = mean, geom = "point", shape = 18, size = 4, color = "red") +
labs(
title = "Piacon eltöltött napok száma ATC csoportonként",
x = "ATC csoport",
y = "Napok száma"
) +
theme_minimal() +
theme(legend.position = "none")
A revíziók számának leíró statisztikáját az 1. táblázat foglalja össze. A revíziók átlagos száma a négy ATC csoportban 14.1 és 17.7 között mozgott. A legmagasabb átlagot a J csoport mutatta (M = 17.7, SD = 14.2), míg a legalacsonyabbat az A csoport (M = 14.1, SD = 9.5).
A piacon töltött idő tekintetében az L csoport gyógyszerei átlagosan rövidebb ideje vannak a piacon (M = 3130 nap, SD = 2455), mint az A (M = 4182 nap, SD = 2524) és J (M = 4300 nap, SD = 2774) csoport készítményei.
A violin plotok (3. és 4. ábra) vizuálisan is szemléltetik a csoportok közötti különbségeket.
# Függő változó: revision_number
## Shapiro-Wilk teszt - revision_number
top4_atc %>%
group_by(atc_name) %>%
summarise(
W = round(shapiro.test(revision_number)$statistic, 3),
p_value = format_p(shapiro.test(revision_number)$p.value))
## # A tibble: 4 × 3
## atc_name W p_value
## <chr> <dbl> <chr>
## 1 "A - Tápcsatorna\nés anyagcsere" 0.959 < 0.001
## 2 "J - Szisztémás\nfertőzésellenes szerek" 0.922 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.87 < 0.001
## 4 "N - Központi\nidegrendszer" 0.919 < 0.001
## Lilliefors teszt - revision_number
top4_atc %>%
group_by(atc_name) %>%
summarise(
D = round(lillie.test(revision_number)$statistic, 3),
p_value = format_p(lillie.test(revision_number)$p.value))
## # A tibble: 4 × 3
## atc_name D p_value
## <chr> <dbl> <chr>
## 1 "A - Tápcsatorna\nés anyagcsere" 0.108 < 0.001
## 2 "J - Szisztémás\nfertőzésellenes szerek" 0.138 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.132 < 0.001
## 4 "N - Központi\nidegrendszer" 0.123 < 0.001
# Függő változó: days_on_market
## Shapiro-Wilk teszt - days_on_market
top4_atc %>%
group_by(atc_name) %>%
summarise(
W = round(shapiro.test(days_on_market)$statistic, 3),
p_value = format_p(shapiro.test(days_on_market)$p.value))
## # A tibble: 4 × 3
## atc_name W p_value
## <chr> <dbl> <chr>
## 1 "A - Tápcsatorna\nés anyagcsere" 0.965 < 0.001
## 2 "J - Szisztémás\nfertőzésellenes szerek" 0.945 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.911 < 0.001
## 4 "N - Központi\nidegrendszer" 0.958 < 0.001
## Lilliefors teszt - days_on_market
top4_atc %>%
group_by(atc_name) %>%
summarise(
D = round(lillie.test(days_on_market)$statistic, 3),
p_value = format_p(lillie.test(days_on_market)$p.value))
## # A tibble: 4 × 3
## atc_name D p_value
## <chr> <dbl> <chr>
## 1 "A - Tápcsatorna\nés anyagcsere" 0.069 0.025
## 2 "J - Szisztémás\nfertőzésellenes szerek" 0.094 < 0.001
## 3 "L - Daganatellenes és\nimmunomoduláns szerek" 0.116 < 0.001
## 4 "N - Központi\nidegrendszer" 0.062 0.091
# Ferdeség és csúcsosság
top4_atc %>%
group_by(atc_name) %>%
summarise(
skewness_revision = round(skewness(revision_number, na.rm = TRUE), 3),
kurtosis_revision = round(kurtosis(revision_number, na.rm = TRUE), 3),
skewness_days = round(skewness(days_on_market, na.rm = TRUE), 3),
kurtosis_days = round(kurtosis(days_on_market, na.rm = TRUE), 3))
## # A tibble: 4 × 5
## atc_name skewness_revision kurtosis_revision skewness_days kurtosis_days
## <chr> <dbl> <dbl> <dbl> <dbl>
## 1 "A - Tápcsato… 0.357 2.19 0.286 2.12
## 2 "J - Szisztém… 0.93 3.34 0.378 2.02
## 3 "L - Daganate… 1.55 6.48 0.912 3.05
## 4 "N - Központi… 1.07 4.00 0.528 2.73
# Vizuális módszerek
## Hisztogram - revision_number
ggplot(top4_atc, aes(x = revision_number, fill = atc_name)) +
geom_histogram(bins = 30, color = "white") +
facet_wrap(~atc_name, scales = "free") +
labs(
title = "Revíziók számának eloszlása ATC csoportonként",
x = "Revíziók száma",
y = "Gyakoriság") +
theme_minimal() +
theme(legend.position = "none")
## Q-Q plot - revision_number
ggplot(top4_atc, aes(sample = revision_number, color = atc_name)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~atc_name, scales = "free") +
labs(
title = "Q-Q plot: Revíziók száma ATC csoportonként",
x = "Elméleti kvantilisek",
y = "Minta kvantilisek") +
theme_minimal() +
theme(legend.position = "none")
## Hisztogram - days_on_market
ggplot(top4_atc, aes(x = days_on_market, fill = atc_name)) +
geom_histogram(bins = 30, color = "white") +
facet_wrap(~atc_name, scales = "free") +
labs(
title = "Piacon eltöltött napok eloszlása ATC csoportonként",
x = "Napok száma",
y = "Gyakoriság") +
theme_minimal() +
theme(legend.position = "none")
## Q-Q plot - days_on_market
ggplot(top4_atc, aes(sample = days_on_market, color = atc_name)) +
stat_qq() +
stat_qq_line() +
facet_wrap(~atc_name, scales = "free") +
labs(
title = "Q-Q plot: Piacon eltöltött napok ATC csoportonként",
x = "Elméleti kvantilisek",
y = "Minta kvantilisek") +
theme_minimal() +
theme(legend.position = "none")
A paraméteres statisztikai eljárások alkalmazhatóságának előfeltételeként a függő változók normalitását Shapiro-Wilk és Lilliefors (Kolmogorov-Smirnov korrekció) tesztekkel vizsgáltam.
A revíziók számának eloszlása a Shapiro-Wilk teszt alapján egyik ATC csoportban sem követte a normális eloszlást (p < .001). A Lilliefors teszt eredményei megerősítették ezt a mintázatot (p < .001). A ferdeségi mutatók pozitív értékeket vettek fel (skewness = 0.36–1.55), ami jobbra ferde eloszlásokra utal. A csúcsossági értékek (kurtosis = 2.19–6.48) az L csoportban kifejezetten leptokurtikus eloszlást jeleztek.
A piacon töltött idő esetében hasonló eredmények lárhatók. A Shapiro-Wilk teszt valamennyi csoportban szignifikáns eltérést mutatott a normalitástól (p < .001). A Lilliefors teszt az N csoport kivételével (p = .091) szintén szignifikáns eredményeket adott.
A Q-Q plotok (5-8. ábra) vizuálisan is alátámasztják a normalitástól való eltérést, különösen az eloszlások szélein megfigyelhető szisztematikus eltérésekkel.
Tekintettel a normalitás sérülésére, a csoportok összehasonlításához nem-paraméteres statisztikai eljárásokat alkalmaztam.
# Függő változó: revison_number
## Kruskal-Wallis teszt - revision_number
top4_atc %>%
kruskal_test(revision_number ~ atc_group) %>%
mutate(p = format_p(p))
## # A tibble: 1 × 6
## .y. n statistic df p method
## <chr> <int> <dbl> <int> <chr> <chr>
## 1 revision_number 1103 9.39 3 0.025 Kruskal-Wallis
## Kruskal-Wallis hatásméret (eta²) - revision_number
top4_atc %>%
kruskal_effsize(revision_number ~ atc_group, ci = TRUE)
## # A tibble: 1 × 7
## .y. n effsize conf.low conf.high method magnitude
## * <chr> <int> <dbl> <dbl> <dbl> <chr> <ord>
## 1 revision_number 1103 0.00581 -0.00073 0.02 eta2[H] small
## Dunn-teszt Holm korrekcióval - revision_number
dunn_revision <- dunn_test(revision_number ~ atc_group, data = top4_atc, p.adjust.method = "holm")
dunn_revision %>%
mutate(p.adj = format_p(p.adj))
## # A tibble: 6 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr> <chr>
## 1 revision_number A J 194 247 1.60 0.110 0.552 ns
## 2 revision_number A L 194 432 -1.05 0.292 0.833 ns
## 3 revision_number A N 194 178 0.0537 0.957 0.957 ns
## 4 revision_number J L 247 432 -3.06 0.00219 0.013 *
## 5 revision_number J N 247 178 -1.50 0.133 0.552 ns
## 6 revision_number L N 432 178 1.09 0.278 0.833 ns
## Cliff's delta függvény 95% CI-vel
calc_cliffs_delta <- function(data, var, group1, group2) {
x <- data %>% filter(atc_group == group1) %>% pull({{var}})
y <- data %>% filter(atc_group == group2) %>% pull({{var}})
result <- effectsize::cliffs_delta(x, y, ci = 0.95)
tibble(
group1 = group1,
group2 = group2,
cliffs_delta = round(result$r, 3),
ci_low = round(result$CI_low, 3),
ci_high = round(result$CI_high, 3)
)
}
## Cliff's delta páronként - revision_number
pairs <- combn(unique(top4_atc$atc_group), 2, simplify = FALSE)
map_df(pairs, ~calc_cliffs_delta(top4_atc, revision_number, .x[1], .x[2]))
## # A tibble: 6 × 5
## group1 group2 cliffs_delta ci_low ci_high
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 L A -0.055 -0.152 0.043
## 2 L J -0.137 -0.224 -0.048
## 3 L N -0.059 -0.158 0.042
## 4 A J -0.094 -0.201 0.014
## 5 A N -0.001 -0.118 0.116
## 6 J N 0.088 -0.023 0.197
# Függő változó: days_on_market
## Kruskal-Wallis teszt - days_on_market
top4_atc %>%
kruskal_test(days_on_market ~ atc_group) %>%
mutate(p = format_p(p))
## # A tibble: 1 × 6
## .y. n statistic df p method
## <chr> <int> <dbl> <int> <chr> <chr>
## 1 days_on_market 1103 48.2 3 < 0.001 Kruskal-Wallis
## Kruskal-Wallis hatásméret (eta²) - days_on_market
top4_atc %>%
kruskal_effsize(days_on_market ~ atc_group, ci = TRUE)
## # A tibble: 1 × 7
## .y. n effsize conf.low conf.high method magnitude
## * <chr> <int> <dbl> <dbl> <dbl> <chr> <ord>
## 1 days_on_market 1103 0.0412 0.02 0.07 eta2[H] small
## Dunn-teszt Holm korrekcióval - days_on_market
dunn_days <- dunn_test(days_on_market ~ atc_group, data = top4_atc, p.adjust.method = "holm")
dunn_days %>%
mutate(p.adj = format_p(p.adj))
## # A tibble: 6 × 9
## .y. group1 group2 n1 n2 statistic p p.adj p.adj.signif
## <chr> <chr> <chr> <int> <int> <dbl> <dbl> <chr> <chr>
## 1 days_on_market A J 197 251 0.105 9.17e-1 1 ns
## 2 days_on_market A L 197 443 -5.15 2.58e-7 < 0.… ****
## 3 days_on_market A N 197 178 -0.551 5.82e-1 1 ns
## 4 days_on_market J L 251 443 -5.71 1.13e-8 < 0.… ****
## 5 days_on_market J N 251 178 -0.683 4.94e-1 1 ns
## 6 days_on_market L N 443 178 4.33 1.50e-5 < 0.… ****
## Cliff's delta páronként - days_on_market
map_df(pairs, ~calc_cliffs_delta(top4_atc, days_on_market, .x[1], .x[2]))
## # A tibble: 6 × 5
## group1 group2 cliffs_delta ci_low ci_high
## <chr> <chr> <dbl> <dbl> <dbl>
## 1 L A -0.256 -0.344 -0.163
## 2 L J -0.254 -0.336 -0.169
## 3 L N -0.23 -0.323 -0.133
## 4 A J -0.013 -0.12 0.095
## 5 A N 0.04 -0.077 0.156
## 6 J N 0.047 -0.064 0.157
A négy ATC csoport közötti különbséget Kruskal-Wallis teszttel vizsgáltam. Az eredmények szignifikáns különbséget mutattak a csoportok között (H(3) = 9.39, p = .025). A hatásméret ugyanakkor kicsi volt (η² = 0.006, 95% CI [-0.001, 0.020]). A 95% CI alsó határa 0 közelébe esik (a becslési eljárás miatt numerikusan enyhén negatív), ami a hatás elhanyagolható mértékével konzisztens. A Dunn-féle post-hoc teszt Holm-korrekcióval egyetlen szignifikáns páronkénti különbséget tárt fel. A J csoport (szisztémás fertőzésellenes szerek) szignifikánsan magasabb revíziószámmal rendelkezett, mint az L csoport (daganatellenes szerek), (padj = .013). A Cliff’s delta hatásméret mutató eredménye alapján kis hatásméret írható le a két csoport között (d = -0.137, 95% CI [-0.224, -0.048]).
A piacon töltött idő tekintetében a Kruskal-Wallis teszt erősen szignifikáns különbséget mutatott a csoportok között (H(3) = 48.2, p < .001). A hatásméret kis-közepes mértékű volt (η² = 0.041, 95% CI [0.020, 0.070]). A Dunn-féle post-hoc teszt három szignifikáns páronkénti különbséget tárt fel (padj < .001). Az L csoport gyógyszerei szignifikánsan rövidebb ideje vannak a piacon, mint az A csoport (d = -0.256), a J csoport (d = -0.254) és az N csoport (d = -0.230) készítményei. Ez a mintázat összhangban van az onkológiai területen tapasztalható intenzívebb kutatás-fejlesztési aktivitással és a gyorsabb gyógyszercserélődéssel.
# A therapeutic_area változó tördelése
drugs_clean <- drugs %>%
separate_rows(therapeutic_area, sep = "\\s*;\\s*") %>%
mutate(therapeutic_area = trimws(therapeutic_area))
# N ATC csoport kiszűrése (központi idegrendszer)
drugs_N <- drugs_clean %>%
filter(!is.na(atc_code)) %>%
mutate(atc_group = str_extract(atc_code, "^[A-Z]")) %>%
filter(atc_group == "N") %>%
filter(!is.na(therapeutic_area))
# Összefoglaló statisztika
drugs_N %>%
summarise(
n_drugs = n_distinct(medicine_name),
n_areas = n_distinct(therapeutic_area),
n_substances = n_distinct(active_substance),
n_companies = n_distinct(marketing_authorisation_holder_company_name))
## # A tibble: 1 × 4
## n_drugs n_areas n_substances n_companies
## <int> <int> <int> <int>
## 1 185 70 101 111
A hálózatelemzés a központi idegrendszeri (N) ATC csoportra fókuszált, amely a pszichiátriai és neurológiai gyógyszereket foglalja magában. A vizsgálat célja a terápiás területek és hatóanyagok közötti kapcsolatok feltérképezése volt.
Az N csoportban összesen 185 egyedi gyógyszer, 70 terápiás terület, 101 hatóanyag és 111 forgalmazó cég volt azonosítható. A therapeutic_area változó pontosvesszővel elválasztott többszörös értékeket tartalmazott, amelyeket külön sorokra bontottam a hálózatépítéséhez.
# Csak a gyakoribb terápiás területek és hatóanyagok (min. 3 előfordulás)
top_areas_N <- drugs_N %>%
count(therapeutic_area) %>%
filter(n >= 3) %>%
pull(therapeutic_area)
top_substances_N <- drugs_N %>%
count(active_substance) %>%
filter(n >= 3) %>%
pull(active_substance)
# Szűrt adatok
drugs_N_filtered <- drugs_N %>%
filter(therapeutic_area %in% top_areas_N | active_substance %in% top_substances_N)
# Élsúlyok hozzáadása (hány gyógyszer kapcsolja össze őket)
edges_weighted <- drugs_N_filtered %>%
count(therapeutic_area, active_substance, name = "weight") %>%
rename(from = therapeutic_area, to = active_substance)
# ATC alcsoport kinyerése (N01, N02, stb.)
atc_subgroups <- drugs_N_filtered %>%
select(therapeutic_area, active_substance, atc_code) %>%
mutate(atc_subgroup = str_extract(atc_code, "^N\\d{2}")) %>%
distinct(active_substance, atc_subgroup) %>%
group_by(active_substance) %>%
slice(1) %>%
ungroup()
# Csúcsok létrehozása típussal és ATC alcsoporttal
nodes_areas_f <- drugs_N_filtered %>%
distinct(therapeutic_area) %>%
rename(name = therapeutic_area) %>%
mutate(type = "Terápiás terület", atc_subgroup = NA_character_)
nodes_substances_f <- drugs_N_filtered %>%
distinct(active_substance) %>%
rename(name = active_substance) %>%
mutate(type = "Hatóanyag") %>%
left_join(atc_subgroups, by = c("name" = "active_substance"))
nodes_f <- bind_rows(nodes_areas_f, nodes_substances_f)
# Hálózat létrehozása súlyozott élekkel
network_f <- tbl_graph(nodes = nodes_f, edges = edges_weighted, directed = FALSE)
# Fokszám és klaszter hozzáadása
network_f <- network_f %>%
activate(nodes) %>%
mutate(
degree = centrality_degree(),
cluster = as.factor(group_louvain()))
# Hálózat vizualizációja - Interaktív ábra (jobb átláthatóság)
layout_coords <- create_layout(network_f, layout = "fr")
layout_coords$x <- layout_coords$x * 2
layout_coords$y <- layout_coords$y * 2
p <- ggraph(layout_coords) +
geom_edge_link(aes(width = weight), alpha = 0.3, color = "gray50") +
geom_point_interactive(aes(
x = x, y = y,
size = degree,
color = ifelse(type == "Terápiás terület", "Terápiás terület", atc_subgroup),
alpha = degree,
tooltip = paste0(name, "\nKapcsolatok: ", degree),
data_id = name
)) +
scale_size_continuous(range = c(2, 12)) +
scale_alpha_continuous(range = c(0.5, 1)) +
scale_edge_width_continuous(range = c(0.3, 3)) +
scale_color_manual(
values = c(
"Terápiás terület" = "steelblue",
"N01" = "#E41A1C",
"N02" = "#FF7F00",
"N03" = "#4DAF4A",
"N04" = "#984EA3",
"N05" = "#F781BF",
"N06" = "#A65628",
"N07" = "#999999"
),
na.value = "coral"
) +
labs(
title = "Terápiás területek és hatóanyagok hálózata (Interaktív)",
subtitle = "N (központi idegrendszer) ATC csoport - Kérem, vigye az egeret a pontokra a részletekért",
color = "Típus / ATC alcsoport",
size = "Kapcsolatok száma",
edge_width = "Közös gyógyszerek"
) +
theme_graph() +
theme(legend.position = "right") +
guides(alpha = "none")
girafe(
ggobj = p,
width_svg = 14,
height_svg = 10,
options = list(
opts_zoom(max = 5),
opts_hover(css = "fill:red;stroke:black;stroke-width:2;"),
opts_tooltip(css = "background-color:white;padding:8px;border-radius:5px;border:1px solid gray;font-size:12px;")))
# Alapvető hálózati statisztikák
network_stats <- tibble(
Mutató = c(
"Csúcsok száma (nodes)",
"Élek száma (edges)",
"Sűrűség (density)",
"Átlagos fokszám (mean degree)",
"Átmérő (diameter)",
"Átlagos úthossz (mean distance)",
"Tranzitivitás (clustering coefficient)"
),
Érték = c(
vcount(network_f),
ecount(network_f),
round(edge_density(network_f), 3),
round(mean(degree(network_f)), 2),
diameter(network_f),
round(mean_distance(network_f), 2),
round(transitivity(network_f), 3)))
network_stats
## # A tibble: 7 × 2
## Mutató Érték
## <chr> <dbl>
## 1 Csúcsok száma (nodes) 110
## 2 Élek száma (edges) 108
## 3 Sűrűség (density) 0.018
## 4 Átlagos fokszám (mean degree) 1.96
## 5 Átmérő (diameter) 32
## 6 Átlagos úthossz (mean distance) 6.61
## 7 Tranzitivitás (clustering coefficient) 0
# Centralitás mutatók hozzáadása a hálózathoz
network_f <- network_f %>%
activate(nodes) %>%
mutate(
degree = centrality_degree(),
betweenness = round(centrality_betweenness(), 2),
closeness = round(centrality_closeness(), 4),
eigenvector = round(centrality_eigen(), 3))
# Top 20 csúcs különböző centralitás mutatók alapján
network_f %>%
activate(nodes) %>%
as_tibble() %>%
arrange(desc(degree)) %>%
select(name, type, degree, betweenness, closeness, eigenvector) %>%
head(20)
## # A tibble: 20 × 6
## name type degree betweenness closeness eigenvector
## <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 Parkinson Disease Terá… 14 167. 0.0345 1
## 2 amifampridine phosphate Ható… 13 150 0.0312 0
## 3 Schizophrenia Terá… 11 69.5 0.0625 0
## 4 Epilepsy Terá… 10 176. 0.02 0
## 5 Migraine Disorders Terá… 6 15 0.167 0
## 6 Alzheimer Disease Terá… 5 54.3 0.0204 0.234
## 7 Opioid-Related Disorders Terá… 5 10 0.2 0
## 8 Pain Terá… 5 51.5 0.0172 0
## 9 duloxetine Ható… 5 93 0.0145 0
## 10 Bipolar Disorder Terá… 4 13.5 0.0333 0
## 11 Epilepsies, Partial Terá… 4 41.5 0.0125 0
## 12 sodium oxybate Ható… 4 12 0.125 0
## 13 Sleep Initiation and Maintena… Terá… 3 5 0.2 0
## 14 Depressive Disorder, Major Terá… 3 41 0.0116 0
## 15 Cancer Terá… 3 78.5 0.0263 0
## 16 Neuralgia Terá… 3 63 0.0167 0
## 17 Narcolepsy Terá… 3 9 0.111 0
## 18 Pain, Postoperative Terá… 3 3 0.333 0
## 19 rivastigmine Ható… 3 37 0.0263 0.365
## 20 pregabalin Ható… 3 118 0.0189 0
A kétoldalú (bipartite) hálózat 110 csúcspontot (terápiás területek és hatóanyagok) és 108 élt (kapcsolatot) tartalmaz. A hálózat ritkának tekinthető (sűrűség = 0.018), ami azt jelenti, hogy az összes lehetséges kapcsolat mindössze 1.8%-a realizálódik. Az átlagos fokszám 1.96, ami arra utal, hogy egy csúcspont átlagosan közel két másik csúcsponttal áll kapcsolatban.
A hálózat átmérője 32, az átlagos úthossz pedig 6.61, ami viszonylag szétszórt, kevéssé integrált struktúrára utal. A tranzitivitás (klaszterezési együttható) értéke 0, ami a bipartite hálózatok sajátossága, ahol azonos típusú csúcspontok nem kapcsolódnak közvetlenül egymáshoz.
A fokszám-centralitás alapján a legkiemelkedőbb terápiás területek a Parkinson-kór (fokszám = 14), a Szkizofrénia (fokszám = 11) és az Epilepszia (fokszám = 10). Ezek a betegségterületek számos különböző hatóanyaggal állnak kapcsolatban, ami a terápiás lehetőségek sokféleségét tükrözi.
A közöttiség-centralitás (betweenness) tekintetében az Epilepszia (176.0) és a Parkinson-kór (167.0) kiemelkedik, ami arra utal, hogy ezek a terápiás területek híd szerepet töltenek be a hálózatban, összekapcsolva más, egyébként elszigetelt klasztereket.
A hatóanyagok közül az amifampridine-foszfát rendelkezik a legmagasabb fokszámmal (13), amely a Lambert-Eaton myastheniás szindróma és más neuromuszkuláris betegségek kezelésében használatos. A duloxetin (fokszám = 5) és a pregabalin (fokszám = 3) szintén kiemelkedő pozíciót foglal el, ami széles terápiás alkalmazhatóságukat tükrözi (fájdalom, depresszió, szorongás, neuropátia).
# Regressziós modell előkészítése
model_df <- top4_atc %>%
transmute(
revision_number,
days_on_market,
atc_group,
generic,
orphan_medicine,
accelerated_assessment
) %>%
filter(
!is.na(revision_number),
!is.na(days_on_market),
!is.na(atc_group)
) %>%
mutate(
atc_group = factor(atc_group),
generic = factor(generic),
orphan_medicine = factor(orphan_medicine),
accelerated_assessment = factor(accelerated_assessment),
days_on_market_z = as.numeric(scale(days_on_market)))
skimr::skim(model_df)
| Name | model_df |
| Number of rows | 1026 |
| Number of columns | 7 |
| _______________________ | |
| Column type frequency: | |
| factor | 4 |
| numeric | 3 |
| ________________________ | |
| Group variables | None |
Variable type: factor
| skim_variable | n_missing | complete_rate | ordered | n_unique | top_counts |
|---|---|---|---|---|---|
| atc_group | 0 | 1 | FALSE | 4 | L: 418, J: 245, A: 190, N: 173 |
| generic | 0 | 1 | FALSE | 2 | FAL: 837, TRU: 189 |
| orphan_medicine | 0 | 1 | FALSE | 2 | FAL: 931, TRU: 95 |
| accelerated_assessment | 0 | 1 | FALSE | 2 | FAL: 996, TRU: 30 |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| revision_number | 0 | 1 | 15.62 | 12.45 | 0.00 | 6.00 | 13.00 | 22.00 | 89.00 | ▇▃▁▁▁ |
| days_on_market | 0 | 1 | 3886.63 | 2542.55 | 124.00 | 1816.50 | 3462.50 | 5488.25 | 9968.00 | ▇▇▆▂▃ |
| days_on_market_z | 0 | 1 | 0.00 | 1.00 | -1.48 | -0.81 | -0.17 | 0.63 | 2.39 | ▇▇▆▂▃ |
A revíziók számát befolyásoló tényezők azonosítására többváltozós lineáris regressziós modelleket építettem. A függő változó a revíziók száma volt, a prediktorok pedig az ATC csoport, a piacon töltött idő (z-standardizált), a generikus státusz, az orphan gyógyszer státusz és a gyorsított értékelés.
Az elemzés 1026 megfigyelésen alapult a hiányzó értékek kizárása után. A piacon töltött idő változót z-standardizáltam az együtthatók értelmezhetősége érdekében.
Két modellt hasonlítottam össze: hagyományos legkisebb négyzetek (OLS) regressziót és robusztus regressziót (MM-becslő), amely kevésbé érzékeny a kiugró értékekre.
# Regressziós modellek létrehozása - két modell: lm vs. lmrob (robusztus)
mod_lm <- lm(
revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment,
data = model_df)
mod_rob <- lmrob(
revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment,
data = model_df)
# Regressziós modellek eredményei
summary(mod_lm)
##
## Call:
## lm(formula = revision_number ~ atc_group + days_on_market_z +
## generic + orphan_medicine + accelerated_assessment, data = model_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -33.208 -4.405 -0.668 4.441 62.413
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 13.7840 0.7468 18.458 < 2e-16 ***
## atc_groupJ 3.1040 0.9543 3.253 0.00118 **
## atc_groupL 3.6263 0.8707 4.165 3.38e-05 ***
## atc_groupN 2.2845 1.0407 2.195 0.02838 *
## days_on_market_z 7.2030 0.3300 21.825 < 2e-16 ***
## genericTRUE -4.1782 0.8287 -5.042 5.46e-07 ***
## orphan_medicineTRUE -1.0943 1.1333 -0.966 0.33449
## accelerated_assessmentTRUE 3.6924 1.8587 1.986 0.04725 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.828 on 1018 degrees of freedom
## Multiple R-squared: 0.3808, Adjusted R-squared: 0.3765
## F-statistic: 89.42 on 7 and 1018 DF, p-value: < 2.2e-16
summary(mod_rob)
##
## Call:
## lmrob(formula = revision_number ~ atc_group + days_on_market_z + generic +
## orphan_medicine + accelerated_assessment, data = model_df)
## \--> method = "MM"
## Residuals:
## Min 1Q Median 3Q Max
## -37.543 -3.935 -0.299 4.045 60.442
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.2598 0.4997 28.536 < 2e-16 ***
## atc_groupJ 3.6720 0.8605 4.267 2.16e-05 ***
## atc_groupL 2.9631 0.6137 4.828 1.59e-06 ***
## atc_groupN 1.5387 0.6807 2.261 0.02399 *
## days_on_market_z 8.8971 0.4217 21.099 < 2e-16 ***
## genericTRUE -3.5214 0.4465 -7.886 7.98e-15 ***
## orphan_medicineTRUE -0.1621 0.5224 -0.310 0.75643
## accelerated_assessmentTRUE 3.6881 1.1926 3.092 0.00204 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Robust residual standard error: 6.065
## Multiple R-squared: 0.6258, Adjusted R-squared: 0.6232
## Convergence in 20 IRWLS iterations
##
## Robustness weights:
## 28 observations c(33,45,238,257,293,302,303,305,411,501,976,990,991,992,993,994,995,1014,1015,1018,1019,1020,1021,1022,1023,1024,1025,1026)
## are outliers with |weight| = 0 ( < 9.7e-05);
## 103 weights are ~= 1. The remaining 895 ones are summarized as
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.001371 0.828100 0.953900 0.850500 0.986900 0.999000
## Algorithmic parameters:
## tuning.chi bb tuning.psi refine.tol
## 1.548e+00 5.000e-01 4.685e+00 1.000e-07
## rel.tol scale.tol solve.tol zero.tol
## 1.000e-07 1.000e-10 1.000e-07 1.000e-10
## eps.outlier eps.x warn.limit.reject warn.limit.meanrw
## 9.747e-05 4.351e-12 5.000e-01 5.000e-01
## nResample max.it best.r.s k.fast.s k.max
## 500 50 2 1 200
## maxit.scale trace.lev mts compute.rd fast.s.large.n
## 200 0 1000 0 2000
## psi subsampling cov
## "bisquare" "nonsingular" ".vcov.avar1"
## compute.outlier.stats
## "SM"
## seed : int(0)
# Együtthatók összehasonlítása
coef_comparison <- bind_rows(
tidy(mod_lm, conf.int = TRUE) %>% mutate(model = "OLS"),
tidy(mod_rob, conf.int = TRUE) %>% mutate(model = "Robust"))
coef_comparison %>%
select(model, term, estimate, std.error, conf.low, conf.high, p.value) %>%
mutate(
estimate = round(estimate, 3),
std.error = round(std.error, 3),
conf.low = round(conf.low, 3),
conf.high = round(conf.high, 3),
p.value = format_p(p.value))
## # A tibble: 16 × 7
## model term estimate std.error conf.low conf.high p.value
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr>
## 1 OLS (Intercept) 13.8 0.747 12.3 15.2 < 0.001
## 2 OLS atc_groupJ 3.10 0.954 1.23 4.98 0.001
## 3 OLS atc_groupL 3.63 0.871 1.92 5.34 < 0.001
## 4 OLS atc_groupN 2.28 1.04 0.242 4.33 0.028
## 5 OLS days_on_market_z 7.20 0.33 6.56 7.85 < 0.001
## 6 OLS genericTRUE -4.18 0.829 -5.80 -2.55 < 0.001
## 7 OLS orphan_medicineTRUE -1.09 1.13 -3.32 1.13 0.334
## 8 OLS accelerated_assessmentT… 3.69 1.86 0.045 7.34 0.047
## 9 Robust (Intercept) 14.3 0.5 13.3 15.2 < 0.001
## 10 Robust atc_groupJ 3.67 0.86 1.98 5.36 < 0.001
## 11 Robust atc_groupL 2.96 0.614 1.76 4.17 < 0.001
## 12 Robust atc_groupN 1.54 0.681 0.205 2.87 0.024
## 13 Robust days_on_market_z 8.90 0.422 8.07 9.72 < 0.001
## 14 Robust genericTRUE -3.52 0.447 -4.40 -2.65 < 0.001
## 15 Robust orphan_medicineTRUE -0.162 0.522 -1.19 0.862 0.756
## 16 Robust accelerated_assessmentT… 3.69 1.19 1.35 6.03 0.002
# Együtthatók vizualizációja
ggplot(coef_comparison %>% filter(term != "(Intercept)"),
aes(x = estimate, y = term, color = model)) +
geom_point(position = position_dodge(width = 0.5), size = 3) +
geom_errorbarh(aes(xmin = conf.low, xmax = conf.high),
position = position_dodge(width = 0.5), height = 0.2) +
geom_vline(xintercept = 0, linetype = "dashed", color = "gray50") +
labs(
title = "Regressziós együtthatók összehasonlítása",
subtitle = "OLS vs. Robusztus regresszió",
x = "Becsült együttható (95% CI)",
y = "Prediktor",
color = "Modell"
) +
theme_minimal()
Az OLS modell szignifikánsan illeszkedett az adatokra, F(7, 1018) = 89.42, p < .001, és a revíziószám varianciájának 38.1%-át magyarázta (R² = .381, R²adj = .377).
A piacon töltött idő volt a legerősebb prediktor: minden egységnyi növekedés a z-standardizált piacon töltött időben 7.20 egységgel növelte a revíziók számát, B = 7.20, 95% CI [6.56, 7.85], p < .001. Ez azt jelenti, hogy a régebben engedélyezett gyógyszerek több revízión estek át.
Az ATC csoport szignifikáns prediktornak bizonyult. A referenciakategóriához (A csoport) képest a J csoport (B = 3.10, p = .001), az L csoport (B = 3.63, p < .001) és az N csoport (B = 2.28, p = .028) egyaránt szignifikánsan magasabb revíziószámmal rendelkezett.
A generikus státusz szignifikáns negatív összefüggést mutatott a revíziók számával (B = -4.18, 95% CI [-5.80, -2.55], p < .001), ami arra utal, hogy a generikus gyógyszerek átlagosan 4.18-cal kevesebb revízión esnek át, mint az originális készítmények.
A gyorsított értékelés szignifikáns pozitív prediktora volt a revíziószámnak (B = 3.69, 95% CI [0.05, 7.34], p = .047). Az orphan gyógyszer státusz nem mutatott szignifikáns összefüggést (B = -1.09, p = .334).
A robusztus regresszió (MM-becslő) 28 megfigyelést azonosított kiugró értékként. A robusztus modell lényegesen magasabb magyarázóerővel rendelkezett (R² = .626, R²adj = .623), ami arra utal, hogy a kiugró értékek jelentősen torzították az OLS becsléseket.
A prediktorok mintázata hasonló maradt, de néhány fontos eltérés mutatkozott. A piacon töltött idő hatása erősebb volt a robusztus modellben (B = 8.90, 95% CI [8.07, 9.72], p < .001). Illetve a gyorsított értékelés hatása is robusztusabbá vált (B = 3.69, p = .002).
A keresztvalidáció (80-20% train-test split) alapján a robusztus modell kissé jobb prediktív teljesítményt mutatott: RMSE = 11.5 vs. 11.7, MAE = 7.32 vs. 7.68, R²pred = .376 vs. .356.
A Cook-féle távolság alapján 51 megfigyelés (5.0%) minősült potenciális kiugró értéknek (Cook’s d > 4/n). A residuumok vizsgálata heteroszkedaszticitásra utalt, ami a robusztus módszerek alkalmazását indokolttá tette.
# Modell illeszkedési mutatók
compare_performance(mod_lm, mod_rob, metrics = c("R2", "R2_adjusted", "RMSE", "AIC", "BIC"))
## # Comparison of Model Performance Indices
##
## Name | Model | R2 | RMSE | AIC (weights) | BIC (weights)
## -----------------------------------------------------------------
## mod_lm | lm | 0.381 | 9.789 | 7610.9 (>.999) | 7655.3 (>.999)
## mod_rob | lmrob | 0.626 | 9.956 | (>.999) | (>.999)
# Regressziós modellek keresztvalidációja
set.seed(123)
idx <- sample(seq_len(nrow(model_df)), size = floor(0.8 * nrow(model_df)))
train <- model_df[idx, ]
test <- model_df[-idx, ]
m1 <- lm(
revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment,
data = train)
m2 <- lmrob(
revision_number ~ atc_group + days_on_market_z + generic + orphan_medicine + accelerated_assessment,
data = train)
pred1 <- predict(m1, newdata = test)
pred2 <- predict(m2, newdata = test)
# Predikciós hibák
rmse <- function(y, yhat) sqrt(mean((y - yhat)^2))
mae <- function(y, yhat) mean(abs(y - yhat))
r2_pred <- function(y, yhat) 1 - sum((y - yhat)^2) / sum((y - mean(y))^2)
prediction_comparison <- tibble(
Modell = c("OLS", "Robust (lmrob)"),
RMSE = round(c(rmse(test$revision_number, pred1), rmse(test$revision_number, pred2)), 3),
MAE = round(c(mae(test$revision_number, pred1), mae(test$revision_number, pred2)), 3),
R2_pred = round(c(r2_pred(test$revision_number, pred1), r2_pred(test$revision_number, pred2)), 3))
prediction_comparison
## # A tibble: 2 × 4
## Modell RMSE MAE R2_pred
## <chr> <dbl> <dbl> <dbl>
## 1 OLS 11.7 7.68 0.356
## 2 Robust (lmrob) 11.5 7.32 0.376
# Regressziós modellek diagnosztikája (ábra)
check_model(mod_lm)
check_model(mod_rob)
# Cook's distance az OLS modellhez
model_df$cooks_d <- cooks.distance(mod_lm)
model_df$fitted <- fitted(mod_lm)
model_df$residuals <- residuals(mod_lm)
# Kiugró értékek azonosítása (Cook's d > 4/n)
n <- nrow(model_df)
outliers <- model_df %>%
filter(cooks_d > 4/n)
cat("Kiugró értékek száma:", nrow(outliers), "\n")
## Kiugró értékek száma: 51
cat("Kiugró értékek aránya:", round(nrow(outliers)/n * 100, 1), "%\n")
## Kiugró értékek aránya: 5 %
# Vizualizáció
ggplot(model_df, aes(x = fitted, y = residuals)) +
geom_point(aes(size = cooks_d, color = cooks_d > 4/n), alpha = 0.6) +
geom_hline(yintercept = 0, linetype = "dashed") +
scale_color_manual(values = c("FALSE" = "gray50", "TRUE" = "red")) +
labs(
title = "Residuumok vs. Illesztett értékek",
subtitle = "Piros pontok: kiugró értékek (Cook's d > 4/n)",
x = "Illesztett értékek",
y = "Residuumok",
size = "Cook's d",
color = "Kiugró"
) +
theme_minimal()
# Összefoglaló táblázat (modelsummary)
modelsummary(
list("OLS" = mod_lm, "Robust" = mod_rob),
estimate = "{estimate} ({std.error})",
statistic = "p.value",
stars = TRUE,
title = "Regressziós modellek összehasonlítása",
coef_rename = c(
"atc_groupJ" = "ATC: J",
"atc_groupL" = "ATC: L",
"atc_groupN" = "ATC: N",
"days_on_market_z" = "Piacon töltött idő (z-score)",
"genericTRUE" = "Generikus",
"orphan_medicineTRUE" = "Orphan gyógyszer",
"accelerated_assessmentTRUE" = "Gyorsított eljárás"
)
)
| OLS | Robust | |
|---|---|---|
| + p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001 | ||
| (Intercept) | 13.784 (0.747) | 14.260 (0.500) |
| (<0.001) | (<0.001) | |
| ATC: J | 3.104 (0.954) | 3.672 (0.860) |
| (0.001) | (<0.001) | |
| ATC: L | 3.626 (0.871) | 2.963 (0.614) |
| (<0.001) | (<0.001) | |
| ATC: N | 2.285 (1.041) | 1.539 (0.681) |
| (0.028) | (0.024) | |
| Piacon töltött idő (z-score) | 7.203 (0.330) | 8.897 (0.422) |
| (<0.001) | (<0.001) | |
| Generikus | -4.178 (0.829) | -3.521 (0.447) |
| (<0.001) | (<0.001) | |
| Orphan gyógyszer | -1.094 (1.133) | -0.162 (0.522) |
| (0.334) | (0.756) | |
| Gyorsított eljárás | 3.692 (1.859) | 3.688 (1.193) |
| (0.047) | (0.002) | |
| Num.Obs. | 1026 | 1026 |
| R2 | 0.381 | 0.626 |
| R2 Adj. | 0.377 | 0.623 |
| AIC | 7610.9 | |
| BIC | 7655.3 | |
| Log.Lik. | -3796.426 | |
| F | 89.422 | |
| RMSE | 9.79 | 9.96 |
Jelen elemzés az Európai Gyógyszerügynökség gyógyszer-engedélyezési adatbázisát vizsgálta négy fő kutatási kérdés mentén.
Engedélyezési trendek. A gyógyszer-engedélyezések száma jelentős fluktuációt mutatott 1995 és 2023 között, kiugró növekedéssel a COVID-19 pandémia időszakában. Az L csoport (daganatellenes szerek) dominanciája különösen az elmúlt évtizedben vált hangsúlyossá.
ATC csoportok közötti különbségek. Szignifikáns különbségek mutatkoztak a revíziók számában (H(3) = 9.39, p = .025) és a piacon töltött időben (H(3) = 48.2, p < .001) a négy leggyakoribb ATC csoport között. Az L csoport gyógyszerei szignifikánsan rövidebb ideje vannak a piacon, ami feltételezhetően az onkológiai területen tapasztalható gyorsabb innovációval magyarázható.
Revíziószámot befolyásoló tényezők. A regressziós elemzés alapján a piacon töltött idő, az ATC csoport, a generikus státusz és a gyorsított értékelés szignifikáns prediktorai a revíziók számának. A robusztus regresszió magasabb magyarázóerőt mutatott (R² = .63 vs. .38), ami a kiugró értékek jelentős torzító hatására utal.
Hálózati kapcsolatok. A központi idegrendszeri gyógyszerek hálózatelemzése feltárta a terápiás területek és hatóanyagok közötti komplex kapcsolatrendszert. A Parkinson-kór, a Szkizofrénia és az Epilepszia központi pozíciót foglal el a hálózatban, ami ezeknek a betegségterületeknek a terápiás diverzitását tükrözi.
Az elemzés korlátai közé tartozik, hogy az adatbázis csak az EMA által centralizált eljárásban engedélyezett gyógyszereket tartalmazza, így az eredmények nem általánosíthatók a nemzeti szinten engedélyezett készítményekre. A revíziók száma mint kimeneti változó korlátozottan tükrözi a gyógyszerek tényleges biztonságossági és hatékonysági profilját. A hálózatelemzés keresztmetszeti jellegű, és nem tárja fel az időbeli változásokat a kapcsolatrendszerben.
További kutatások vizsgálhatnák a revíziók tartalmát és típusát, a forgalomból történő kivonások prediktorait, valamint az engedélyezési folyamat időtartamának meghatározó tényezőit.